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SECTION  I 


INTRODUCTION 


The  treatment  of  soil-structure  interaction  is  of  considerable  importance  in 
analyses  of  the  integrity  of  structures  in  ground-shock  environments.  There  is 
currently  one  basic  approach  to  the  nonlinear  treatment  of  the  problem:  finite 
element  or  finite  difference  methods.  Finite  element  or  finite  difference  methods 
can  model  the  problem  to  almost  any  accuracy  desired  but  the  large  number  of  equa¬ 
tions  required  generally  precludes  efficient  computation.  The  number  of  equations 
can  be  reduced  through  special  elements  representing  quiet  boundaries  or  specific 
radiation  properties,  but  modeling  with  these  special  elements  requires  consider¬ 
able  skill  and  may  necessitate  a  complete  finite  element  model  for  validation.  An 
approach  to  achieve  a  method  more  versatile  than  special  elements  and  more  economi¬ 
cal  than  finite  element  or  finite  difference  for  the  treatment  of  these  problems 
may  be  to  consider  a  boundary  element  formulation.  Such  an  approach  is  pursued  in 
this  study:  an  analytical  approximation  of  the  soil-structure  interaction  requir¬ 
ing  boundary  elements  is  combined  with  the  modeling  capabilities  of  the  finite 
element  (FE)  method  while  avoiding  the  burden  of  many  elements  in  the  soil. 

This  report  examines  a  boundary-element  (BE)  treatment  of  the  surrounding  non¬ 
linear  soil  that  is  an  extension  of  a  BE  method  for  a  linear  soil  [1],  In  the  lin¬ 
ear  problem,  the  structure  was  modeled  with  an  available  FE  code,  and  the  soil- 
structure  interaction  was  reduced  to  a  surface  relationship  through  the  use  of  a 
doubly  asymptotic  approximation  (DAA)  [2],  which  required  the  application  of  BE 
techniques  [3].  The  linear  study  focused  on  the  two-dimensional  plane-strain  re¬ 
sponse  of  structures  surrounded  by  an  infinite  elastic  medium;  the  approach  was 
shown  to  produce  results  of  acceptable  accuracy. 

In  the  nonlinear  soil  problem,  the  linear  structure  is  still  modeled  with  an 
available  FE  code.  The  soil-structure  interaction  is  reduced  to  a  surface  relation¬ 
ship  for  the  linear  portion  of  the  soil  behavior  and  a  volume  relationship  for  the 
nonlinear  portion  through  the  use  of  a  modified  DAA,  which  still  involves  the  appli¬ 
cation  of  BE  techniques.  The  nonlinear  problem  requires  a  considerable  increase  in 
the  complexity  of  the  analysis.  One,  volume  information  is  needed  to  evaluate  the 
volume  relationship;  hence,  additional  data  must  be  generated  and  stored.  Two,  the 
DAA  must  be  modified  to  provide  a  wave  propagation  model  for  soil  response  based  on 
an  Ij-plasticity  model.  And  three,  the  wave-decoupling  that  occurs  from  assuming  the 
incident  and  scattered  waves  are  algebraically  additive  must  be  examined  carefully. 
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To  study  the  nonlinear  problem  effectively,  a  simple  problem  that  retains  the 
physics  but  eliminates  much  of  the  generality  was  developed.  The  simple  problem  is 
the  axisymmetric  response  of  an  internally  loaded,  plane  strain  circular  infinite  shell 
surrounded  by  an  infinite  2-D  nonlinear  soil.  This  model  eliminates  the  wave¬ 
decoupling,  since  there  is  no  incident  wave,  but  includes  the  first  two  aspects  of 
the  nonlinear  problem.  A  simple  study  of  the  wave-decoupling  effect,  not  considered 
in  this  report,  is  given  in  Appendix  A. 

The  report  first  addresses  the  outline  of  the  theory:  the  linear  theory  is 
reviewed,  the  governing  nonlinear  equations  are  given,  and  alternate  approaches  are 
considered.  Then  the  computational  strategy  is  presented,  followed  by  numerical 
results  and  conclusions.  Figures  and  tables  are  grouped  at  the  end  of  the  section 
in  which  they  are  first  mentioned. 


SECTION  II 


TECHNICAL  DISCUSSION 

The  first  part  of  this  section  describes  the  theoretical  development  of  doubly 
asymptotic,  boundary- integral  analysis  techniques  for  nonlinear  media.  First,  the 
elastic  formulation  is  reviewed  [1,  2].  Then  the  extension  to  inelastic  soil 
response  is  described.  The  introduction  of  the  doublv  asymptotic  approximation  leads 
to  two  possible  formulations  for  the  treatment  of  inelastic  effects;  both  of  these 
are  considered.  In  addition,  two  alternate  formulations  are  presented. 

Treatment  of  an  inelastic  material  by  the  boundary  element  method  has 
been  described  in  considerable  detail  by  Swedlow  and  Cruse  [3)  and  by  Mendelson  [4], 
Mukerjee  [5],  however,  has  shown  that  these  descriptions  are  incorrect  for  two-dimen¬ 
sional  plane-strain  response.  The  results  of  these  previous  developments  are  used  in 
this  report  without  detailed  derivation;  however,  an  outline  of  the  correct  deriva¬ 
tion  is  given  in  Appendix  B. 

The  second  part  of  this  section  describes  the  computational  strategy  used  to 
solve  the  equations  governing  inelastic  soil-structure  interaction.  Empha¬ 

sis  is  placed  on  the  methods  used  to  determine  the  volume  information,  displacements, 
strains  and  nonlinear  effects.  The  details  of  the  linear  portion  of  the  computations 
have  been  treated  previously  [ 1 } . 

2.1  REVIEW  OF  THE  LINEAR  DAA  FORMULATION 

The  matrix  FE  equation  of  motion  for  a  linear-elastic  structure,  embedded  in  a 
surrounding  linear-elastic  medium  and  excited  by  known  forces  applied  to  the  struc¬ 
ture  and  by  an  incident  wave  propagating  through  the  medium,  is 

Mq+Kq=f+f+fc  (1) 

~s  -  ~s  3  -s-I-S  ' 

where  M  and  K  are  the  mass  and  stiffness  matrices  for  the  structure,  q  is  the 
structural  displacement  vector,  fg  is  a  vector  of  known  forces  applied  to  the  struc¬ 
ture,  fj  and  fg  are  surface-force  vectors  associated  with  the  incident  and  scattered 
waves,  respectively,  and  a  dot  denotes  temporal  differentiation.  The  DAA  is  now 
introduced  to  evaluate  the  scattered  force  vector  fg  [2].  This  approximation,  which 
is  a  surface  interaction  approximation  that  replaces  the  infinite  volume  of  an  ex¬ 
ternal  surface  of  the  structure,  is  expressed  as  [1] 

f„  =  p  DT  GT  A  C  G  uc  +  DT  K  u c  (2) 

-S  ~  ~  ~  ~m  ~  -  S  ~  ~m  -S 
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where  p  is  the  mass  density  of  the  medium,  £  and  G  are  coordinate  transformation 
matrices,  A  is  a  diagonal  surface  element  area  matrix,  is  a  diagonal  sound-speed 
matrix  for  the  medium,  K  is  a  full  surface-stiffness  matrix  for  the  medium,  u„  is 

~iT\  — b 

the  computational  scattered-displacement  vector  for  the  surface  elements,  and  a 
superscript  T  denotes  matrix  transposition.  The  doubly  asymptotic  nature  of  the 
DAA  is  apparent  from  (2),  i.e.,  at  low  frequencies,  u^  is  small  relative  to  u  ,  and 
(2)  reduces  to  a  static  stiffness  relation;  at  high  frequencies,  the  reverse  is 
true,  and  (2)  reduces  to  a  radiation  damping  relation. 


With  structure-medium  surface  compatibility  and  elastic  field  superposition 
requiring  that  Dq  =  u  +  u_ ,  (1)  and  (2)  may  be  combined  to  obtain 

“1  *"b 


t  T 

M  q+pD'AC  G  D  q  +  (K  +D  K  D) 
~s  -  ~  ~  ~m  ~  ~  -  ~s  ~  ~m  ~ 


q  = 


f  +  fT  +  p  DT  GT  A  C  G  ut  +  DT  K  u t 
-s  -I  ~m  ~  -I  ~  -in  -I 


(3) 


In  the  development  of  the  nonlinear  problem  that  follows,  only  known  forces  applied 
to  the  structure  are  considered,  so  f ^  =  u^  =  0.  For  this  situation,  (3)  reduces  to 


M  q+pDTGTAC  G  D  q  +  (K  +  DT  K  D)q  =  f 


■'S  -*  ~  ^  ~ 

2.2  NONLINEAR  DAA  FORMULATION 


-  s 


(4) 


In  (4),  the  elements  of  C  are  the  dilitational  and  shear  sound  speeds  for  the 

rn 

medium  [1].  One  effect  of  inelastic  behavior  is  to  alter  these  sound  speeds;  at 
present,  however,  this  alteration  is  neglec  ,;d  on  the  assumption  that  the  dominant 
effect  of  inelastic  behavior  of  the  medium  manifests  itself  in  a  displacement- 
dependent  form.  Also  note  that  the  assumptions  (u  =  u  +  u  and  f  =  fT  +  fQ)  alter 
the  functional  dependence  of  the  incident  and  scattered  waves,  an  important  consid¬ 
eration  for  the  nonlinear  problem.  A  simple  study  of  these  assumptions  are  pre¬ 
sented  in  Appendix  A,  but  the  effects  of  the  assumptions  are  not  considered  further 
in  this  report. 

In  order  to  include  the  displacement-dependent  effects  of  material  inelasticity, 
it  is  necessary  to  consider  the  two-dimensional  boundary-integral  equation  [3,4,3] 


Xuk(P) 


+ 


+ 


(P,Q)u*(Q)dL(Q)  = 


L 

J C*mk(p  ,Q)e*m(p)dA(p)  - 
A 


(P,Q)t*(Q)dL(Q) 
(p ,Q)b£ (p)dA(p) 


A 


(5) 
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where  P  is  a  point  either  on  the  structure-medium  interface  L  or  in  the  surrounding 

region  A,  Q is  a  point  on  L,  p  is  a  point  in  A,  A  =  ^  if  P  is  on  L  and  A  =  1  if  P  is 
k.  It 

in  A,  u  and  t  are  components  of  the  displacement  and  traction  vectors  pertaining 

kJ.  k? 

to  the  kth  Cartesian  direction,  respectively,  T  and  U  are  components  of  second- 

order  tensors,  which  constitute  Green's  functions,  E  is  a  component  of  a  third- 

•  ?.in 

order  tensor,  which  also  constitutes  a  Green's  function,  e  is  a  component  of  the 

Z  ^ 

plastic-strain  tensor  and  b  is  a  component  of  a  body  force  vector  in  A,  Through  the  div¬ 
ision  of  L  into  a  series  of  boundary  elements  (BE),  and  A  into  an  array  of  finite 
quadrature  elements  (QE)  ,  (5)  may  be  expressed  in  matrix  notation,  for  u  and  t  on  L,  as 

Su  =  Ft  +  Be  -  Eb  (6) 


in  which  the  2x2  elements  of 

S,  F  and  E 

are  given  by 

II 

o*  *n 

in 

^iAt  + 

f  Tk?  £?  dL. 

J  *3  3  3 

L  . 

3 

^  ■ 

f  ukt  n. 

J  IJ  3 

L  . 

1 

dL  . 

3 

■ 

/  "o'  *i 

dA  . 

3 

(7) 

j 

and  the  2x3  elements  of  B  are  given  by 


^  Amk 

ij 


/_  A  m  k  A  m  , 
l .  .  n .  dA  . 
ij  J  J 


(8) 


A  . 
J 


In  (6)  Cp  is  arranged  in  groups  as  (ep 


11  12 


b22). 


In  (7), 


<5^  are  Kronecker  deltas,  i  and  j  are  BE  indices  for  L  in  the  first  two  of 


-p  >  ep2)j  and  -  as  (b11 

6 .  .  and 

^  ‘  l  ? 

(7),  j  is  a  QE  index  for  A  in  the  third  of  (7),  £  and  q  are  assumed  BE  shape- 

Z 

functions,  is  an  assumed  QE  shape-function,  L.  is  the  length  of  the  jth  boundary 

ks,  k£  3 

element,  and  the  kernels  T. .  and  U..  are  given  in  Appendix  B.  In  (8),  i  is  a  BE 


ij 


index  for  L,  j  is  a  QE  index  for  A, 
area  of  the  jth  QE,  and  the  kernel 


Am 


is  an  assumed  QE  shape-function,  A.  is  the 


Amk 


is  given  in  Appendix  B. 


Am 


£  p  £ 

In  the  present  implementation,  the  shape  functions  £  ^ ,  n^ ,  ifr ^ ,  and  q 

taken  as  unity  over  each  element,  so  that  each  element  is  described  by  a  single, 

kA 

centrally  located,  nodal  point.  Hence,  constitutes  the  displacement  in  the  kth 
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direction  at  nodal  point  i  due  to  a  unit  point  load  in  the  &th  direction  at  a  point 

j^nik. 

on  element  j.  The  associated  stresses  at  i  are  given  by  o-h  =  e  where  the  e, 

k£  Omk  1,  1j  *  k  * 

are  the  Cartesian  base  vectors.  Incidentally,  T. .  =  I..  n*5,  where  the  n  are  the 

il  13  1  1 

direction  cosines  for  the  surface  normal  at  i.  The  numerical  integration  method 
used  to  evaluate  the  first  two  of  (7)  is  described  in  [1],  and  the  method  used  for 
the  last  of  (7)  and  for  (8)  is  given  in  Appendix  C. 


In  (6)  F  is  a  nonsingular  matrix;  hence  (6)  may  be  rewritten  as 

tc  =  K  u  -  F~ 1  B  e  +  F~ 1  E  b  (9) 

-S  ~m  -S~~-p~~- 

where  K  =  F  *  S  is  the  linear  surface-stiffness  matrix  for  the  medium  [see  (2)] 
~m  ~  ~ 

and  the  subscript  S  has  been  added  to  indicate  that  the  scattered  wave  is  modeled. 
Equations  (2)  and  (9)  are  now  combined  to  give 


fc  =  p  DT  GT  A  C  G  u.  +  DT  K  uc  -  DT  F~ 1  Be  +  DT  F_1  E  b  (10) 

which  represents  a  possible  DAA  for  an  inelastic  medium.  In  the  development  that 
follows,  the  body  force  b  is  not  considered  until  the  alternate  formulations  are 
presented;  see  Section  2.3.  In  2.3.1  the  body  force  is  viewed  as  the  acceleration 
of  the  medium  and  the  linear  damping  term,  the  first  term  in  (10),  is  replaced  by  the 
body  force  term.  The  requirements  and  assumptions  used  to  derive  (3)  and  (4)  are 
applied  again  to  combine  (1)  and  (10)  to  obtain  the  governing  nonlinear  equation  of 
motion 


M  q  + 


T  T 

p  D  G  A  C 
~  ~  ~  ~ra 


G  D  q  +  (K  +  DT  K  D)  q 
~  ~  —  ~s  ~  ~m  1 


+  Q 


(11) 


where  Q  =  £*  F  *  B  is  a  psuedo-force  vector  that  accounts  for  medium  inelas¬ 
ticity.  As  the  DAA  is  a  surface  approximation  [1]  and  the  evaluation  of 
Q  requires  volume  information  to  calculate  £  ,  an  additional  model  must  be  con¬ 
structed  to  generate  the  volume  information.  Two  possible  volume-information  models 
are  now  examined. 


2.2.1  Quasi-Static  DAA 

Inasmuch  as  the  psuedo-force  vector  Q  may  be  viewed  as  a  correction  to  the 

T 

linear  surface-stiffness  force-vector  D  ^  £  q,  a  possible  volume-information 
model  might  derive  from  the  static  solution  for  the  medium  based  on  the  current 
soil-structure  tractions  and  displacements.  In  [3,4,5]  integral-equation  formulas 
are  given  for  medium  strains  due  to  boundary  tractions  and  displacements;  these  may 
be  cast  into  matrix  notation  for  easy  computational  implementation.  Unfortunately, 
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these  strain  formulas  did  not  prove  to  be  computationally  efficient.  The  kernels 

1^.  £ 

were  poorly  behaved  because  they  are  based  on  derivatives  of  the  U _  kernel,  which 
itself  exhibits  marginal  numerical  behavior.  This  difficulty  was  overcome  through 
the  use  of  the  boundary-integral  equation  for  displacements  due  to  boundary  trac¬ 
tions  and  displacements,  i.e.,  (5)  with  A  =  1.  In  matrix  notation  (5)  becomes,  for 
A  -  1, 


u 

-m 


=  -  Su  +  Ft  +  Be 

- - - P 


(12) 


where  u^  is  the  displacement  vector  in  the  medium,  u  and  t  are  the  boundary-displace¬ 
ment  and  boundary- traction  vectors,  the  2x2  elements  of  S  and  F  are  given  by 


,k£  _  f 

'lj  =  / 


L. 

J 


■  / 


Tk!"  dl. 

1 J  J  J 


ki  Z 

U  .  .  n  .  dL . 

3-3  J  J 


(13) 


and  the  2x3  elements  of  jB  are  given  by 


|  imk 

ij 


/ 

Ai 


„£mk  im  , 

I  .  .  n  .  d  A. 

3-3  3  3 


(H) 


Most  of  the  symbols  in  (12),  (13),  and  (14)  have  been  defined  previously;  here, 
however,  i  is  a  QE  index  for  A  and  j  is  a  BE  index  for  L  in  (13),  and  i  and  j  are 
QE  indices  for  A  in  (14). 

Equation  (12)  efficiently  provides  the  displacements  in  the  medium;  from  these 
displacements,  the  strains  are  quickly  computed  from  finite-difference  expressions. 
With  the  strains  available,  the  computations  to  determine  from  the  inelastic 
medium  behavior  proceed  in  a  straightforward  manner;  the  details  are  given  in  Sec¬ 
tion  2.4. 

A  difficulty  was  encountered  in  using  the  quasi-static  DAA  formulation  for  one 
of  the  two  inelas tic  soil  models  considered .  To  explain  this  difficulty,  it  is  neces¬ 
sary  to  discuss  briefly  the  two  inelastic  soil  models. 

The  first  model  considered,  the  mechanical  sublayer  model  (6] ,  is  a  basic  J ^  (second 
deviatoric  stress-tensor  invariant)  flow  rule  theory  commonly  used  for  the  elasto- 
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plastic  behavior  of  metals.  This  model  was  chosen  because  it  is  simple  and  because 
a  computational  subroutine  was  readily  available  at  the  laboratory.  It  was  recog¬ 
nized  that  this  model  would  not  properly  represent  soil  behavior  but  it  would  pro¬ 
vide  a  first  step  in  the  verification  of  the  nonlinear  DAA  formulation.  As  shown 
in  Section  3.0,  the  quasi-static  DAA  with  this  plasticity  model  produced  acceptable 
results. 

The  second  model  considered  is  a  cap  model  [7],  herein  called  the  Cap  model  (also 
see  [8]  for  details  of  the  Cap  model  subroutine).  This  inelastic  soil  model  was  de¬ 
veloped  to  model  granular  soils  under  dynamic  loading  conditions  and  is  based  on  a 
yield  (failure)  surface  that  is  a  function  of  1^  (first  stress-tensor  invariant)  and 

Although  this  model  is  regarded  as  an  excellent  model  for  dynamically  loaded 
soils,  it  appears  to  be  incompatible  with  the  static  soil  behavior  assumed  for  the 
quasi-static  DAA  formulation. 

To  understand  the  difficulty  encountered  with  the  Cap  model,  the  static  solu¬ 
tion  of  the  2-D  plane  strain  cylindrical  cavity  of  radius  a  for  internal  loading  in 
an  infinite  elastic  medium  is  considered.  For  an  internal  axisymmetric  applied  dis¬ 
placement,  u  ,  or  a  pressure  giving  rise  to  u  ,  the  displacement  and  stresses  in  the 
a  a 

elastic  medium  are 

u  a 

a  _  „  _ 


and 

2 u  aG 
a 

a  =  - ~ — 

r  2 

r 


e  ’ 


(16) 


where  u  and  u  are  the  radial  and  tangential  displacements,  r  is  the  radial  coor- 
r  0 

dinate  measure,  a,  and  a  are  the  stresses  in  cylindrical  coordinates  and  G  is 

T  0  Z 

the  elastic  shear  modulus.  With  this  stress  field  (16),  the  expression  for  1^  be¬ 
comes 


1 1  H  ar  +  0e  +  °z  =  0  (17) 

Hence,  for  this  simple  axisymmetric  static  problem,  Ij  is  always  zero,  and  the  quasi¬ 
static  DAA  formulation  produces  a  stress  field  with  no  Ij  component.  In  the  actual 
dynamic  problem,  inertial  forces  in  the  medium  give  rise  to  a  significant  Ij  compon¬ 
ent,  as  well  as  a  ^  component,  so  that  the  stress  state  at  the  failure  surface  is 
quite  different  from  that  predicted  by  the  quasi-static  DAA.  Hence,  the  quasi-static 
DAA  formulation  does  not  produce  responses  that  are  comparable  to  a  finite-element 
dynamic  response  calculation. 
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In  summary,  although  the  result  for  the  quasi-static  DAA  with  a  plasticity 
model  were  encouraging,  the  results  with  the  Cap  model  were  not.  Therefore  a  DAA 
formulation  that  more  realistically  reflects  the  dynamic  nature  of  the  problem  was 
explored. 

2.2.2  Quasi-Dynamic  DAA 

As  the  DAA  was  known  to  produce  good  results  for  an  elastic  problem,  the  physi¬ 
cal  model  that  corresponds  to  the  elastic  DAA  was  sought.  This  development  can  most 
easily  be  described  in  terms  of  the  following  problem:  the  plane-strain  excitation 
of  an  infinite  medium  surrounding  an  infinite,  circular  cylindrical  cavity  by  means 
of  transient  pressurization  of  the  cavity. 

For  a  linear  medium,  the  DAA  says  that  radial  stress  is  given  by 

ar(r,t)  =  -  P  c*  ur(r,t)-~  ur(r,t)  (18) 

where  the  DAA  surface  (i.e.,  the  surface  defined  by  the  soil-structure  interface) 

has  been  positioned  at  an  arbitrary  radial  position  r.  In  this  equation,  p,  c  and 

<t> 

U  are  the  density,  dilatational  velocity  and  shear  modulus  for  the  medium,  respec¬ 
tively,  u^  is  radial  displacement,  and  the  dot  denotes  temporal  differentiation. 

But,  from  elasticity  theory, 

°r  (r  ,  t )  =  P  if  -  ~  ur(r,t) 

Ur(r,t)  =  97  ^ ( r » c  >  (19) 

where  <f>  is  the  dilatational  displacement  potential.  Hence  (18)  and  (19)  yield,  for 
the  "internal  forcing"  problem  considered  here, 

4»(r,t)  =  *(a,t  -  — )  (20) 

C<P 

where  a  is  the  radius  of  the  cavity.  Now  this  result  constitutes  a  plane-wave  treat¬ 
ment  of  the  radiated  wave.  Hence  the  DAA  not  only  provides,  through  (18),  a  radiated- 
wave  stress-displacement  relation  at  r  =  a,  but  also  provides  the  displacement  field 
in  the  medium,  viz., 

ur  ( r  ,  t )  =  u  (a,t  -  —■)  (21) 

4> 

For  an  elasto-plastic  medium,  a  surface  relation  such  as  (18)  is  not  sufficient 
to  determine  surface  response.  In  this  case,  "volume  information"  is  also  needed. 
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Such  information  is  provided  by  application  of  the  method  of  characteristics  to  the 
radiated  wave,  which  leads  to  the  behavior  illustrated  in  Figure  1,  The  character¬ 
istics,  which  constitute  loci  of  constant  displacement,  define  the  displacement 
field  in  the  medium  at  any  time  of  interest.  Straight  characteristics  indicate 
linear  wave  propagation,  while  curved  characteristics  indicate  nonlinear  wave  pro¬ 
pagation. 


In  accordance  with  Figure  1,  nonlinear  DAA  computations  proceed  as  follows: 


1.  Surface  response  is  determined  at  t  =  At  based  on  linear-elastic  medium 
behavior. 

2.  The  surface  displacement  ur(a,0)  is  prescribed  at  r  -  a  =  c  At  and  average 

$ 

strains  are  calculated  for  the  region  0  <  r  -  a  <  c  it. 

<P 

3.  The  state  of  the  medium  at  r  =  At  is  found  to  be  linear-elastic. 

4.  Surface  response  is  determined  at  t  =  2At  based  on  linear-elastic  medium 
behavior. 


5.  u  (a, 0)  and  u  (a, At)  a>-e  prescribed  at  r  -  a  =  2c  At  and  r  -  a  =  c  At,  re- 

*  L  <P  (f> 

spectively;  average  strains  are  calculated  for  the  region  0  <  r  -  a  <  2c  At. 

6.  The  state  of  the  medium  at  t  =  2At  is  found  to  be  linear-elastic  for 

r  -  a  >  c  At  and  inelastic  for  r  -  a  <  c  At. 

9  <P 

7.  A  local  propagation  veloicty  c^  =  c(r=a)  is  calculated  and  surface  response 
is  determined  at  t  =  3 At  based  on  inelastic  medium  behavior. 


8. 


u  (a,0) ,  u  (a, At)  and  u  (a,2At)  are  prescribed  at  r  -  a  =  3c  At,  r  -  a  = 

L  t  T  m 

and  r  -  a  =  c^At,  respectively;  average  strains  are  calculated  for  the 

region  0  <  r  -  a  £  3c, At. 

9 


2c  At 
9 


9.  The  state  of  the  medium  at  t  =  3At  is  found  to  be  linear-elastic  for 

r  -  a  z  2c, At  and  inelastic  for  r  -  a  <  2c  At. 

9  $ 


10.  Local  propagation  velocities  c^  =  c(r=a)  and  c^  =  cCr^a+c^At)  are  calculated 
and  surface  response  is  determined  at  t  =  4At  based  on  inelastic  medium  be¬ 
havior. 


Note  that  this  procedure,  which  pertains  to  a  volume-wave  (VW)  model,  defines  a  new 
characteristics  grid  at  each  time  step.  For  computational  purposes,  two  single¬ 
dimension  arrays  are  required  to  store  the  displacement  and  its  radial  location;  for 
a  nonaxisymraetric  problem  these  would  be  two-dimensional  arrays.  Also  in  the  cal¬ 
culations,  characteristics-grid  displacements  and  slopes  (first-derivatives)  are 
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1 


K 


interpolated  to  yield  nodal  values  for  a  fixed  spatial  grid,  from  which  medium 
strains  are  calculated.  This  is  required  for  efficient  computation  of  the  volume 
integral  for  nonlinear  DAA  analysis;  otherwise  the  boundary  matrix  j5  [see  02)] 
would  have  to  be  recomputed  at  each  time  step.  Note  the  volume  matrices  S,  F  and 

15  in  (12)  are  not  required,  so  considerably  less  storage  is  required. 

For  the  volume-wave  model  just  discussed,  radial  displacement  does  not  attenu¬ 
ate  with  distance  [see  (21)].  However,  during  the  study  of  an  alternate  formulation, 

as  discussed  in  Section  2.3,  an  attenuation  of  — i.  was  found  to  represent  more  accu- 

/r 

rately  the  exact  response  (as  might  be  expected).  Hence  this  attenuation  was 
actually  used  to  generate  the  characteristics-grid  displacements  for  the  VW  model. 

The  quasi-dynamic  DAA  was  partially  successful.  The  no-I^  problem  encountered 
with  the  quasi-static  DAA  was  not  encountered  with  the  quasi-dynamic  formulation. 

The  characteristics-grid  displacements,  described  above,  agree  with  FE  solutions  for 
very  early  times.  At  later  times,  however,  agreement  fades,  because  the  tensile- 

cutoff  feature  in  the  Cap  model  produces  displacement  responses  in  the 

medium  that  are  not  properly  treated  by  the  simple  characteristics  method  developed 
for  the  quasi-dynamic  DAA.  As  partial  compensation,  the  DAA  characteristics-grid  was 
modified  to  treat  tensile  cutoff  (a  loss  of  soil  strength  in  tension)  by  not  allowing 
a  displacement  to  propagate  into  or  out  of  a  region  of  tensile  cutoff.  The  dis¬ 
placement  in  the  tensile  cutoff  region  was  interpolated  from  the  two  points  adjacent 
to  the  region.  This  means  that  back  reflections  from  the  inner  surface  or  a  region 
in  cutoff  are  ignored.  Such  reflections  would  be  difficult  to  include  in  the  cal¬ 
culations  and  would  border  on  actually  solving  for  the  dynamic  displacement  field  in 
the  medium.  The  goal  of  this  study  has  been  to  avoid  costly  calculations;  the  in¬ 
clusion  of  reflections  would  be  costly. 

2.3  ALTERNATE  FORMULATIONS 

Two  alternate  formulations  of  the  dynamic  soil-structure  interaction  problem  and 
a  modification  of  the  DAA  approach  are  considered  here.  The  first  formulation  treats 
the  radiation  damping  in  terms  of  a  body-force  field  associated  with  the  acceleration 
field  in  the  medium,  which  is  obtained  from  the  volume-wave  model  discussed  in  Section 
2.2.2.  The  second  formulation  is  based  on  the  use  of  "infinite  elements"  [10]  in  lieu 
of  boundary  elmeents.  Finally,  the  DAA  modification  addresses  tensile  cutoff  in  the 
Cap  plasticity  model. 
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2.3.1  Body  Force  Formulation 


This  formulation  is  best  explained  by  examining  the  governing  equation  for  the 
medium  and  then  reconsidering  the  arguments  put  forward  to  develop  the  DM  [1,2]. 

For  two-dimensional  plane  strain,  including  plastic  flow,  the  Navier  equation,  which 
governs  the  behavior  of  the  medium,  may  be  written  as 


V 


+ 


1 

1-2V 


(u 


j  >j 


+ 


— ^ — (r  P  ) 

l-2v  Uk,k'i 


+  p'u. 


i>  j  >k  =  1,2,  (22) 

2 

where  v  is  the  Laplacian  operator,  u  is  the  displacement  field,  v  is  Poisson's 
ratio,  G  is  the  elastic  shear  modulus,  ep  is  the  plastic  strain  field,  b  is  the  body 
force  field,  p  is  the  mass  density  and  a  superimposed  dot  indicates  a  time  deriva¬ 
tive;  see  Appendix  B  for  a  derivation  of  (22).  Terms  appearing  in  the  DM  relation 
(10)  are  readily  associated  with  terms  in  (22).  First,  the  stiffness-force  term  in¬ 
volving  is  associated  with  the  left  side  of  (22);  second,  the  plastic  term  Cp 
is  associated  with  the  Ep  terms  in  (22);  finally,  the  damping  force  term  involving 
is  associated  with  the  inertia  term  pu^  in  (22)  (see,  e.g.,  [9]). 

Rigorously  speaking,  the  inertia  term  is  costly  to  compute,  requiring  a  finite- 
element  mode]  of  the  medium  for  practical  problems.  However,  the  volume-wave  model 
described  in  Subsection  2.2.2  directly  provides  the  displacement  history  of  the  medium, 
so  that  the  acceleration  iiL  may  easily  be  computed  from  a  three-point  backward  dif¬ 
ference  formula.  The  static  boundary  integral  equation  (5)  does  not  admit  an  accel¬ 
eration  term,  but  inertial  forces  may  easily  be  included  by  considering  them  as  body 
forces.  Hence,  from  (22),  the  body-force  computational  vector  in  (9)  is  given  by 

b  =  -  pGu  (23) 

so  that  the  surface  relationship  is  given  by  (10)  with  C  =  0  and  b  is  given  by  (23). 

~'Tn  - 

Hence,  for  f  =  f^  =  0,  (1)  becomes 

M  q  +  (K  +  DT  K  D) q  =  f0  +  Q  -  DT  F_1  E  b.  (241 

~s  -  ~s  ~  ~m  -S  2  ~  ~  ~  ~  v 


Solutions  obtained  from  (24)  were  quite  satisfactory  for  elastic  medium  response, 
but  n  satisfactory  for  inelastic  response.  To  obtain  converged  solutions,  the  med¬ 
ium  quadrature  grid  must  be  uniform  and  fixed- increment  time  integration  is  required. 
Even  then,  the  medium  solution  jumps  about  too  much  to  obtain  smooth  inelastic  re¬ 
sponse. 
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2.3.2  Infinite  Elements 

In  this  section,  some  recent  developments  in  the  application  of  finite-element 
methods  to  infinite  media  are  discussed  [10].  Although  an  alternate  formulation 
based  on  infinite  elements  was  not  pursued  in  this  study,  the  formulation  merits  con¬ 
sideration,  as  it  is  attractive  from  two  viewpoints  [11].  First,  the  method  is  based 
on  variational  principles  that  provide  symmetric  matrices  without  resort  to  ati  hoc 
symmetry  procedures.  Second,  the  extensive  software  developed  for  finite-element 
methods  is  directly  applicable  to  these  infinite  elements.  The  second  item  is  prob¬ 
ably  the  more  important,  as  a  considerable  portion  of  this  study  effort  was  devoted 
to  the  development  of  special  software  for  boundary-integral  equation  methods.  To 
date,  infinite-element  concepts  have  been  app’ied  only  to  fluid  media  [10],  so  fur¬ 
ther  development  work  would  be  required  for  application  to  inelastic  solid  media. 


2.3.3  Separation  Model 

One  additional  modeling  concept  was  employed  in  an  attempt  to  account  more  real¬ 
istically  for  tensile  cutoff  phenomena  embodied  in  the  Cap  plasticity  model  f7,8]. 

For  the  axisymmetric  problems  studied  here,  tensile  cutoff  occurs  on  a  circle  speci¬ 
fied  by  a  certain  radial  distance  out  from  the  soil-structure  interface.  As  tensile 
cutoff  represents  the  inability  of  the  soil  to  support  a  tensile  load,  the  soil  be¬ 
yond  the  tensile  cutoff  circle  no  longer  provides • stiffness  and  separation  occurs. 
Hence,  after  tensile  cutoff  occurs,  the  system  of  interest  is  assumed  to  be  composed 
of  the  structure  and  a  surrounding  annulus  of  soil,  i.e.,  the  soil  is  no  longer  in¬ 
finite.  This  separation  model  can  be  partially,  partially  in  that  reflections  from 
the  tensile  cutoff  surface  are  not  considered,  incorporated  in  the  computations. 

The  incorporation  involves  a  modification  to  the  elastic  soil  contribution,  the 
matrix,  and  the  inelastic  contribution,  Q.  The  procedures  to  affect  the  modifications 
are  presented  below,  first  for  the  elastic  modification  and  then  the  inelastic  one. 


where  p^  is  the  internal  pressure,  ur(a)is  the  radial  displacement  at  the  radius  a 
where  the  pressure  p^  is  applied,  G  is  the  soil's  shear  modulus,  v  is  the  soil's 
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Poisson's  ratio,  and  b  is  the  radius  of  the  tensile  cutoff  circle.  As  a  check,  (26) 
reduces  to  (25)  as  b  -»■  °°.  Hence  the  ratio  of  the  two  soil  stiffnesses  is  obtained  as 

/k. 


=  Ka/Ki 


(27) 


and  the  separation  due  to  tensile  cutoff  is  accounted  for  by  scaling  by  <  to  pro¬ 
duce  the  correct  stiffness  for  the  annulus.  The  calculations  for  Q,  the  inelastic 
component,  due  to  e  are  then  performed  for  the  annulus  only.  Note  that  the  same 
model  would  be  obtained  if  the  soil  matrices  in  (12)  and  (24)  were  recomputed  for 
the  proper  annulus,  but  this  would  be  very  expensive.  Unfortunately,  the  numerical 
results.  Section  3.0,  show  that  this  additional  "refinement"  in  the  model  did  not 
greatly  improve  the  results. 

2.4  COMPUTATIONAL  STRATEGY 


The  governing  equations  of  motion  for  the  DAA  formulation  (11)  and  for  the  body 
force  formulation  (24)  are  simply  nonlinear  second-order  equations  that  are  common 
in  structural  dynamics  analysis.  Therefore,  these  equations  are  readily  solved  by 
direct  time  integration  methods  for  nonlinear  structural  dynamics  analysis  [12,13, 
14].  These  methods  fall  into  one  of  two  categories:  an  implicit  method  [12]  and  an 
explicit  method  [13,14];  their  relative  merits  will  be  discussed  in  Section  3.  Both 
require  that  software  be  available  to  generate  the  linear  coefficient  matrices,  the 
nonlinear  terms  and  the  forcing  function.  The  procedures  for  the  linear  terms  and 
the  forcing  function  are  presented  in  [1],  so  the  emphasis  here  is  on  the  nonlinear 
terms. 


First,  a  brief  overview  of  the  steps  required  to  compute  the  nonlinear  terms  is 
presented.  Then  the  techniques  used  to  implement  the  volume-wave  model  and  the  soil 
inelasticity  treatment  are  discussed. 

2.4.1  Overview  of  Nonlinear  Computations 

The  plastic  strain  eP  is  the  key  quantity  that  must  be  computed  to  solve  the 
nonlinear  equations  (11)  and  (24).  The  medium  acceleration  is  also  needed  for  (24), 
but  this  is  relatively  simple  to  compute  and  will  be  dealt  with  later.  The  following 
steps  are  involved  in  the  computation  of  eP: 

1.  At  the  n-th  time  step,  the  solution  vector  q11  is  known;  hence,  from  u°  =  Dqn, 
un  is  known. 

2.  The  vector  un  is  used  to  compute  the  displacements  in  the  medium,  either 
from  (12)  or  in  accordance  with  the  volume-wave  model  described  in  Subsec¬ 
tion  2.2.2.  The  details  of  this  calculation  are  considered  below. 
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3.  Through  application  of  the  proper  differentiation  formulas,  the  total- 
strain  vector  en  is  determined  for  the  centroids  of  the  quadrature  elements 
in  the  medium.  The  details  of  this  calculation  are  also  considered  below. 

4.  The  strain-increment  vector  is  determined  from  Ae11  =  e”  -  e"  \  where  e11  * 
has  been  retained  from  the  previous  time  step. 

5.  The  strain-increment  vector  Ae11  and  the  previous  total-stress  vector  a11  ^ 
are  sufficient  to  compute  the  total-stress  vector  at  the  n-th  step,  a11, 
based  on  the  inelastic  constitutive  relationships  discussed  in  Subsection 
2.2.1. 

6  n 

6.  The  elastic-strain  vector  (e  )  is  computed  as 


( e 
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n 


n 
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where  is  the  array  of  elastic-material  coefficients  that  relates  strains 
to  stress. 

7.  The  plastic-strain  vector  (eP)n  is  obtained  as 


,p.n  _  n  .en 
(e  )  =  e  -  (e  ) 

Figure  2  illustrates  steps  6  and  7. 

8.  The  nonlinear  force  Qn  is  computed  by  matrix-vector  multiplication;  see  the 
definition  following  (11). 


9a.  For  the  explicit  integrator,  the  solution  at  n+1  is  obtained  from  informa- 

n+l  .  , 

txon  at  n,  so  q  is  now  computed. 

9b.  For  the  implicit  integrator,  the  solution  at  n+1  is  obtained  from  an  extra¬ 
polation  based  on  Qn  and  Qn  The  corrected  value  of  Qn+^ ,  which  is  com¬ 
puted  from  the  solution  gn+^is  then  compared  to  Qn+^  -  extrapolated;  if  the 
difference  lies  within  an  acceptable  error  bound,  the  integration  proceeds; 
if  not,  an  iteration  procedure  is  carried  out  until  satisfactory  converg¬ 
ence  is  achieved  within  a  reasonable  number  of  iteration  cycles  or  the  cal¬ 
culation  is  aborted. 


2.4.2  Volume  Information  Computations 

The  displacements  in  the  medium  are  computed  either  from  a  quasi-static  model, 
as  described  in  Subsection  2.2.1,  or  from  a  quasi-dynamic  model,  as  described  in  Sub¬ 
section  2.2.2.  The  reason  for  computing  the  displacements  in  the  medium  is  to  deter¬ 
mine  strains  by  numerical  differentiation.  The  plastic  strains,  assumed  to  be  uniform 
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over  each  quadrature  element,  are  computed  at  fixed  grid  points  at  the  center  of 
these  prescribed  quadrature  elements.  A  sketch  of  some  typical  quadrature  elements 
and  the  plastic  strain  points  is  shown  in  Figure  3. 


For  the  quasi-static  model  the  displacements  in  the  medium  are  computed  from 
(12),  where  the  coefficient  matrices  are  evaluated  such  that  the  displacements  are 
obtained  at  the  points  marked  with  an x in  Figure  3.  These  points  are  chosen  to  be 
half-way  between  the  plastic  strain  points  in  both  the  radial  and  tangential  direc¬ 
tions.  This  location  is  ideal,  as  it  avoids  singular  points  in  the  evaluation  of 
(12)  and  provides  simple  but  accurate  differentiation  formulas  to  compute  the  strains. 
As  an  example,  the  strains  at  strain  point  4  are  computed  as 
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where  the  subscripts  r,6  refer  to  the  polar-coordinate  directions  and  the  superscripts 
refer  to  the  medium  displacement  indices.  Note  that,  for  strain  points  adjacent  to 
the  boundary,  e.g.,  1,3,5,  ...  ,  the  boundary-displacement  value  is  used  as  one  of 
the  medium  displacement  points.  In  this  study  only  calculations  for  axisymmetric 
problems  were  performed,  so  y  and  the  second  term  in  e  vanished.  For  general  geo- 
metries,  the  quadrature  elements  may  be  viewed  as  finite  elements,  the  medic  <  dis¬ 
placements  computed  at  the  proper  nodes,  and  standard  finite-element  methods  used  to 
compute  the  strain  for  the  quadrature  (finite)  element  [15]. 

For  the  quasi-dynamic  model,  the  medium  displacements  are  not  computed  at  fixed 
points;  instead,  grid  points  are  generated  by  the  location  of  a  particular  displace¬ 
ment  on  its  characteristic.  These  grid  points  are  positioned  as 
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where  d  is  the  radial  distance  from  the  boundary  at  time  t  ,  At.  is  the  i-th  time 
n  n  j  J 

step  and  c(d.)  is  the  dilatational  sound  speed  for  the  medium  at  the  d.  location. 

1  J 

The  variation  in  the  medium's  sound  speed  is  determined  from 


+  j  u 


(30) 


where  B,  the  bulk  modulus,  and  p,  the  shear  modulus,  are  determined  from  the  stress 
and  strain  increments  for  the  quadrature  element  under  consideration. 

As  indicated  in  the  displacement  snapshot  shown  in  Figure  4,  the  displacement 
of  the  medium  at  d^  is  computed  as 


u(tn) 


a  +  d 


(31) 


where  u(t  )  is  the  radial  displacement  at  the  boundary  at  time  t  ,  and  a  is  the  rad- 
n  r  n 

ius  of  the  boundary.  This  corresponds  to  cylindrical  spreading  of  the  scattered  wave. 

To  compute  strains,  the  slope  3u^/3r  and  the  displacement  uf  are  required  at  each 
plastic-strain  point.  These  values  are  determined  as  follows.  For  plastic-strain 
point  1 ,  the  average  slope  and  displacement  over  the  quadrature  element  are  computed 
as  (see  Figure  4) 
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where  u(AR.)  is  evaluated  at  AR,  by  linear  interpolation  between  u  ,  and  u 
1  i  n-1 


n-2  * 


For 


plastic-strain  point  2  in  Figure  4,  the  average  slope  is  simply  the  slope  of  the 

straight  line  joining  u^  and  un_2»  an<*  the  average  displacement  at  the  plastic- 

strain  point  is  determined  by  linear  interpolation  between  u  ,  and  u 

n-1  n-2 

After  the  strains  are  computed,  Steps  4-9  in  Subsection  2.4.1  are  straightfor¬ 
ward,  as  these  steps  are  typical  of  any  transient  structural  dynamics  analysis. 
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[12,13,14,15].  Finally,  the  medium-acceleration  values  required  in  (23)  are  provided 
by  numerical  differentiation  based  upon  displacement  values  at  three  consecutive  time 
steps. 

Note  that  displacements  and  strains  are  computed  in  polar  coordinates  [see,  e.g., 
(32)],  while  the  coefficient  matrices  for  the  governing  equations  (11)  and  (24)  are 
computed  in  Cartesian  coordinates,  which  are  the  global  coordinates.  Therefore, 
polar-coordinate  displacements  and  strains  are  transformed  to  Cartesian  coordinates 
before  the  matrix-vector  multiplications  in  (11)  or  (24)  are  performed. 

Because  the  nonlinear  formulations  require  additional  matrices,  only  one-half 
the  problem  is  modeled  to  reduce  storage  requirements,  with  the  X2~axis  in  Figure  3 
chosen  as  a  plane  of  symmetry.  The  symmetry  condition  is  easily  established  for  the 
FE  structure  [15],  but  that  for  the  BE  model  is  not  as  readily  determined  (see  Appen¬ 
dix  D).  Finally,  it  is  important  to  mention  that,  although  calculations  were  only 
performed  for  axisymmetric  problems,  most  of  the  analysis  software  was  constructed 
to  treat  more  general  non-axisymmetric  problems. 
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SECTION  III 


NUMERICAL  RESULTS 

The  numerical  results  presented  in  this  section  illustrate  the  characteristics 
of  the  soil-analysis  formulations  described  in  Section  2.  All  the  results  pertain 
to  the  axisynunetr ic  response  of  either  an  infinite,  circular  cylindrical  cavity  or 
an  infinite,  circular  cylindrical  shell  surrounded  by  an  infinite  soil-medium  to  an 
internal  triangular  pressure  pulse.  The  computational  model  for  the  shell  is  a  fin¬ 
ite-element  one,  as  provided  by  the  REXBAT  Code  f 1 7 ] .  Two  inelastic  models  for  the 
soil  are  considered:  1)  the  mechanical  sublayer  model  [6]  and  2)  the  Cap  model 
[7,8];  see  Subsection  2.2.1.  As  eight  elements  over  half  the  circumference  (the  x~- 
axis  being  the  symmetry  axis)  produce  excellent  results  for  the  elastic  problem, 
eight  elements  are  used  in  all  the  examples.  The  governing  equations  of  motion  arc 
solved  with  the  stand-alone  time-integration  packages  described  in  [12,  Hand  14]. 
These  packages  permit  both  automatic-step  and  fixed-step  integration  using  the  im¬ 
plicit  Park  method  [12]  or  the  explicit  central-difference  method  [13,14]. 

The  following  comparisons  are  presented  in  this  section: 

1.  Responses  produced  by  the  doubly  asympt<  tic  approximation  (DAA)  are  com¬ 
pared  with  exact  responses  for  an  elastic  medium  surrounding  an  infinite 
cylindrical  cavity; 

2.  DAA  and  body-force  responses  are  compared  with  exact  responses  for  an 
elastic  medium  surrounding  an  infinite  cylindrical  shell; 

3.  Quasi-static  and  quasi-dynamic  DAA  responses  are  compared  with  NONSAP  re¬ 
sponses  for  an  inelastic  (^-plasticity  theory)  medium  surrounding  an 
infinite  cylindrical  cavity; 

4.  Quasi-static  DAA  responses  are  compared  with  TRANAL  responses  for  an  inelas¬ 
tic  (Cap  model)  medium  surrounding  an  infinite  cylindrical  shell; 

5.  Quasi-dynamic  DAA  responses  are  compared  with  TRANAL  responses  for  an  in¬ 
elastic  (Cap  model)  medium  surrounding  an  infinite  cylindrical  shell;  and 

6.  Quasi-dynamic  body-force  responses  are  compared  with  TRANAL  responses  for 
an  inelastic  (Cap  model)  medium  surrounding  an  infinite  cylindrical  shell. 

The  first  three  comparisons  demonstrate  the  satisfactory  accuracy  of  the  DAA 
and  body  force  formulations  for  an  elastic  medium,  and  a  ^-plasticity  model  inelas- 
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tic  medium.  The  remaining  comparisons  demonstrate  the  unsatisfactory  accuracy  of 
these  formulations  for  the  Cap  model  inelastic  media. 

3.1  ELASTIC  MEDIUM:  DAA  AND  EXACT  RESULTS 

In  [1],  the  elastic  response  of  an  infinite  cylindrical  shell  embedded  in  an 
elastic  medium  to  an  incident  dilatational  wave  was  computed  with  the  DAA  formula¬ 
tion  and  compared  with  the  corresponding  exact  response.  Satisfactory  agreement 
was  observed.  In  this  study,  internal  triangular  pressure  pulses  are  considered  in 
lieu  of  an  incident  wave.  Hence,  it  is  appropriate  to  compare  results  produc  ed  bv 
the  DAA  formulation  with  their  exact  counterparts  for  the  internal  loading. 

The  problem  considered  is  the  response  of  a  cavity  in  a  2-D  plane  strain  elastic 
medium  with  an  elastic  shear  modulus,  G,  of  200,000  psi,  Poisson's  ratio,  v,  of  0.25 

3 

and  a  weight  density  of  117  lbs/ft  .  The  internal  pressure  pulse  is  an  isosceles 
triangle  of  durations  t,  2t,  At  and  8t,  where  t  is  the  transit  time  of  a  dilatational 
wave  across  the  diameter  of  the  cavity.  The  exact  response  is  obtained  by  the  resid¬ 
ual  potential  method  [2]  and  the  DAA  response  is  obtained  by  discrete  element  methods 
[1].  The  DAA  peak-displacement  errors  are  20.8%,  24.1%,  22.7%  and  15.7%  for  the 
pulse  widths  of  t,  2t,  At  and  8t,  respectively.  These  DAA  results  are  not  as  accur¬ 
ate  as  those  reported  in  [1],  but  they  are  still  considered  acceptable;  the  introduc¬ 
tion  of  a  cylindrical  shell  would  yield  results  accurate  to  within  10%-20%,  as  ob¬ 
served  in  [1],  Note  that  the  most  inaccurate  responses  are  obtained  for  the  inter¬ 
mediate  pulses,  whose  dominant  frequency  components  lie  in  the  intermediate  frequency 
range,  where  the  DAA  is  most  inaccurate. 

3.2  ELASTIC  MEDIUM:  DAA,  BODY-FORCE  AND  EXACT  RESULTS 

This  problem  was  studied  to  determine  the  proper  volume-wave  model  (Subsection 

2.2.2)  and  to  illustrate  the  accuracy  of  the  body-force  formulation  (Subsection  2.3.1). 

The  problem  models  the  response  of  an  infinite  cylindrical  shell  of  rad ius/thickness 

=  50,  where  the  radius  is  lm  (39.3701  in.),  E  =  30  x  10^  psi,  v  =  0.3,  and  p  =  460 
3 

lb/ft  embedded  in  an  elastic  soil  which  has  the  following  properties:  G  =  559  ksi, 

3 

v  =  0.3313,  and  p  =  125  lb/ft  .  The  pressure  pulse  is  a  triangle  with  a  rise  time 
of  0.07288  ms,  a  duration  of  0.7288  msec  (1  transit  time  across  the  diameter)  and  a 
peak  pressure  of  1  ksi. 

Figure  5  shows  exact  [2]  and  DAA  [see  equation  (A)]  displacement  responses. 

Note  that  the  DAA  underpredicts  the  peak  and  is  more  heavily  damped.  Figure  6  illus¬ 
trates  the  displacement  responses  for  the  exact  and  body-force,  (25)  with  Q  =  0,  form¬ 
ulations.  Note  that  the  peaks  are  nearly  identical;  the  only  difference  occurs  at 
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later  times.  The  difference  occurs  because  the  body-force  calculation  is  very 
sensitive  to  numerical  noise  that  accumulates  throughout  the  computations.  To 
achieve  the  results  shown  in  Figure  6  the  characteristics  volume-wave  model  dis¬ 
placement  was  attenuated  by  /a/r ,  where  a  is  the  radius  of  the  soil-structure  inter¬ 
face  and  r  is  the  distance  from  the  origin.  In  addition  the  true  area  of  the  quad¬ 
rature  elements  was  used.  Without  the  attenuation,  a  cylindrical  plane  wave  model, 
the  DAA  damping  results  are  essentially  reproduced. 

These  results  demonstrate  that  the  volume-wave  model  of  Subsection  2.2.2  is 
highly  accurate  for  elastic  medium  response.  The  body-force  formulation  would 
therefore  be  preferred  over  the  DAA  formulation,  except  that  it  has  the  following 
shortcomings.  First,  the  volume  integration  must  enclose  all  of  the  medium  through 
which  the  wave  front  of  the  outgoing  wave  passes.  Fortunately,  however,  large  ele¬ 
ments  can  be  used;  in  this  problem,  for  example,  18  radial  elements  were  used  to 
reach  10  cavity  radii  out  into  the  medium.  Second,  numerical  noise  produced  by  the 
double-differentiation  of  displacement  was  so  severe  that  obtaining  a  reliable  solu¬ 
tion  was  difficult.  The  only  satisfactory  results  were  obtained  with  quadrature 
elements  of  uniform  radial  dimension  and  time  integration  with  a  small  fixed-step. 
The  time  increment  required  to  obtain  a  satisfactory  solution  was  roughly  one-fourth 
that  required  to  obtain  a  satisfactory  solution  with  the  DAA  formulation.  The  ex¬ 
plicit  central-difference  method  was  much  superior  to  the  implicit  Park  method  for 
solving  the  body-force  equations  of  motion. 

The  numerical  noise  problem  discussed  above  may  be  circumvented  by  other  stra¬ 
tegies;  two  examples  follow.  One,  the  boundary  acceleration,  obtained  from  u  =  Dq , 
could  be  propagated  into  the  medium  as  the  boundary  displacements  are  in  the  VW 
model.  Two,  the  field  acceleration  can  be  computed  from  a  first  differentiation  of 
the  stress  field,  i.e.,  from  a  force  balance.  Either  of  these  methods  should  pro¬ 
duce  a  smoother  acceleration  field  for  computational  purposes.  These  methods  were 
not  pursued  in  this  study. 

3.3  J2-THE0RY  MEDIUM:  DAA  AND  NONSAP  RESULTS 

This  problem  involves  the  elastic-plastic  response  of  an  infinite  cylindrical 
cavity  surrounded  by  a  .^-plasticity  theory  inelastic  medium;  see  section  3.1  for 
the  elastic  response.  The  cavity  has  a  radius  of  23.0  inches.  The  material  proper¬ 
ties  are  E  =  500,000  psi,  v  =  0.25,  C  =  200,000  psi,  a  von  Mises  yield  stress  of 

7500  psi,  a  slope  of  50,000  psi  after  yield  (0.1  the  elastic  slope)  and  a  density 

3 

of  117  lb/ft  .  A  J2  flow  rule  plasticity  model  is  used  for  the  DAA  computations. 
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This  is  based  on  a  two  element  mechanical  sublayer  model  [6].  The  excitation  is 
provided  by  an  internal  pressure  pulse  represented  by  an  isosceles  triangle  with  a 
17,800  psi  peak  and  a  duration  of  1.71  msec  (2  transit  times:  2  diameters).  For 
the  cavity,  Mg  =  Kg  =  JO  so  the  governing  equation  (H)  becomes  a  first-order  equa¬ 
tion.  Therefore,  the  implicit  integration  must  be  used  to  solve  this  problem,  as 
the  central-difference  explicit  integrator  is  valid  only  for  second-order  equations. 

The  radial  displacement  response,  shown  in  Figure  7,  was  computed  with  both 
the  quasi-static  DM  (Subsection  2.2.1)  and  the  quasi-dynamic  DM  (Subsection  2.2.2) 
formulations,  for  comparison  with  the  NONSAP  Finite  Element  Code  [18].  For  the 
quasi-static  DM  calculation,  the  centroids  of  the  volume  (area)  quadrature  elements 
are  located  at  28.75,  40.00,  62.50,  107.50,  and  197.50  inches.  For  the  first  quasi¬ 
dynamic  DM  calculation,  the  centroids  of  the  quadrature  elements  are  at  27.50, 

32.50,  37.50,  42.50,  47.50,  52.50,  67.50,  82.50,  97.50,  112.50,  127,50,  142,50, 

192.50,  242.50,  292.50,  342.50,  392.50  and  442.50  inches.  The  second  quasi-dynamic 
DM  calculation  makes  use  of  only  the  first  6  quadrature  elements  characterizing 
the  first  quasi-dynamic  calculation. 

Figure  7  shows  that  all  the  DM  formulations  produce  similar  responses.  This 
indicates  that  the  quasi-static  and  quasi-dynamic  models  are  essentially  equivalent 
for  an  inelastic  J2  material,  and  that  only  a  relatively  small  volume  near  the  soil- 
structure  interface  need  be  considered  to  obtain  a  converged  solution.  Figure  7 
also  shows  that  the  DM  response  peaks  are  approximately  37%  smaller  than  the  NONSAP 
peak,  and  that  the  DM  permanent  displacements  are  approximately  33%  smaller  than  the 
NONSAP  permanent  displacement.  This  compares  with  the  24.1%  underprediction  of  peak 
response  by  the  DM  in  the  elastic  case  discussed  in  Subsection  3.1. 

This  response  problem  is  considered  to  represent  an  overly  severe  test  of  the 
DM  for  a  material.  This  is  because  no  embedded  structure  is  included  to  miti¬ 
gate  inaccuracies  introduced  by  the  DM  (cf.  Subsection  3.1).  In  any  case,  inelas¬ 
tic  soil  behavior  is  not  described  by  J2  plasticity  theory,  so  the  more  realistic 
Cap  model  must  be  introduced,  as  described  in  the  following. 

3.4  CAP-MODEL  MEDIUM:  QUASI-STATIC  DM  AND  TRANAL  RESULTS 

To  provide  check  cases  for  the  quasi-static  DM  formulation  of  Subsection  2.2.1, 
the  problems  shown  in  Table  1  are  considered.  These  problems  were  run  on  the  DM 
code  and  the  TRANAL  code  (19].  The  material  properties  used  in  these  calculations 
are  shown  in  Table  2.  The  Cap  model  data  [20]  was  modified  slightly  to  provide  an 
elastic  response  before  yielding  by  taking  X0“-R(A-C).  This  modification  was  nec- 
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essary  to  move  the  ^-failure  surface  up  the.  ^-axis;  otherwise  yielding  occurs 
immediately  with  the  quasi-static  DAA  formulation. 

For  these  check  case.s  agreement  between  the  DAA  and  TRANAL  responses  was  non¬ 
existent.  For  the  profile  2  pressure  pulse  the  DAA  response  predicts  Pg  =  4800  psi, 

while  TRANAL  predicts  P  =  1500  psi.  For  the  nonlinear  calculations  at  P  =  2P  and 

e  o  e 

5Pe>  the  agreement  is  even  poorer.  In  addition,  nonlinear  calculations  for  the  same 
peak  pressure  are  in  poor  agreement.  This  poor  agreement  is  characterized  by  magni¬ 
tude  differences  of  a  factor  of  3-10  and  the  shape  of  the  displacement  histories  not 
being  similar.  The  basic  reason  for  the  poor  comparison  is  the  fact  that  the  quasi¬ 
static  DAA  stress  field  has  no  I  component  for  the  internal  axisymmetric  load  (see 
Subsection  2.2.1).  For  this  reason,  the  quasi-dynamic  DAA  formulation  of  Subsection 
2.2.2  was  developed. 

3.5  CAP-MODEL  MEDIUM:  QUASI-DYNAMIC  DAA  AND  TRANAL  RESULTS 

These  check  cases,  shown  in  Table  3,  are  very  similar  to  the  internal-excita¬ 
tion  axisymmetric  cases  considered  in  Section  3.4.  Except  a  slightly  different 
pressure  pulse  is  used;  the  reason  for  this  is  discussed  below.  Now  the  peak  pres¬ 
sure  is  chosen  as  7.35  ksi  instead  of  being  based  on  the  pressure  to  cause  yielding 
of  the  soil.  The  material  properties  are  also  very  similar;  see  Table  4.  Now  xo 
is  -0.0441  ksi  as  given  in  [20].  Since  the  quasi-dynamic  DAA  is  based  on  reproduc¬ 
ing,  as  best  as  possible,  the  actual  volume  displacement,  the  modification  to  to 
move  the  initial  yield  surface  is  no  longer  needed.  And  with  this  xq  there  is 
essentially  no  elastic  response  so  defining  the  peak  based  on  Pg  is  no  longer  applic 
able. 

Figures  8  and  9  show  displacement  response  histories  produced  by  the  quasi¬ 
dynamic  DAA  (Q-D  DAA)  formulation  without  "separation",  along  with  corresponding 
TRANAL  response  histories  [21].  Except  at  very  early  times,  agreement  between  the 
DAA  and  TRANAL  histories  is  nonexistent.  Actually,  the  histories  begin  to  diverge 

when  tensile  cutoff  occurs  in  the  TRANAL  calculations.  This  occurs  roughly  at  t/TA  = 

2 

0.5  for  the  short  pulse  and  t/t^  =  2.0  for  the  long  pulse.  The  crude  DAA  separation 
model  of  Subsection  2.5  hardly  improves  the  situation,  as  shown  in  Figures  10  and  11 
Indeed,  the  DAA  results  of  Figure  11  diverge. 

In  all  of  the  DAA  calculations,  the  volume  (area)  integration  encompassed  a 
region  enclosing  the  farthest  circle  reached  by  the  outgoing  wave  front.  For  the 
short  pulse,  eighteen  radial  elements  were  used,  with  the  centroids  at  40.57,  42.97, 
45.37,  47.77,  50.17,  52.57,  62.17,  71.77,  81.37,  90.97,  100.6,  110.2,  148.6,  187.0, 
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225.4,  263.8,  302.2,  and  340.6  inches.  For  the  long  pressure  pulse,  seventeen  ra¬ 
dial  elements  were  used,  with  the  centroids  at  41.07,  44.47,  47.87,  51.27,  54.67, 
58.07,  73.47,  88.87,  104.3,  119.7,  135.1,  150.5,  220.9,  291.3,  361.7,  432.1,  and 
502.5  inches.  The  quasi-dynamic  Cap-model  calculations  were  quite  sensitive  to 
variations  in  the  volume  mesh;  this  is  in  contrast  to  the  relative  insensitivity 
exhibited  in  the  quasi-dynamic  calculations  for  the  Jj  model. 

The  calculations  presented  in  this  subsection  were  performed  with  the  central- 
difference  explicit  time  integrator.  The  explicit  integrator  was  eventually  chosen 
over  the  implicit  one  because  it  was  faster  and  more  reliable  when  confronted  with 
the  highly  irregular  nonlinear  soil  stiffness  forces  associated  with  the  tensile 
cutoff  and  dilatancy  aspects  of  the  Cap  model. 

3.6  CAP-MODEL  MEDIUM:  BODY-FORCE  RESULTS 

The  last  attempt  to  obtain  reasonable  correlation  with  the  TRANAL  results  in¬ 
volved  the  body-force  formulation  discussed  in  Subsection  2.3.1.  The  problem  solved 
pertains  to  the  short  pressure  pulse  described  in  Table  3.  Body-force  displacement 
response  is  compared  with  TRANAL  response  in  Figures  12  and  13.  The  correlation  be¬ 
tween  body-force  and  TRANAL  results  is  seen  to  be  poor.  A  comparison  of  Figures  8 
and  10  with  Figures  12  and  13,  respectively,  shows  that  the  body-force  and  quasi- 
dynamic-DAA  formulations  produce  similar  responses,  indicating  consistency  between 
the  two  formulations. 

To  obtain  stable  body-force  solutions,  a  small  fixed  time  step  was  used  in  the 
explicit  integrations  and  the  18  volume  elements  were  all  20  inches  in  radial  dimen¬ 
sion.  The  "no  separation"  calculation  involved  a  constant  elastic  sound  speed  in 
the  soil,  while  the  "separation"  calculation  involved  a  variable  sound  speed.  The 
variable  sound  speed  option  produces  the  slowly  growing  instability  shown  in  Figure 
13,  as  discussed  in  Subsection  2.3.1. 


Figure  5. 


Elastic  response:  DAA  and  exact 
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Figure  6.  Elastic  response:  body  force  damping  and  exact 
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Figure  8.  Elastic  plastic  Cap  model  response,  short  pulse 
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Figure  12.  Elastic  plastic  Cap  model  response,  short  pulse,  and  body  force  damping 
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Figure  13.  Elastic  plastic  Cap  model  response,  short  pulse,  body  force  damping 
and  with  "separation". 


Table  1.  DAA-SSI  Check  Runs:  1 


Medium: 


NTS  Tuff  (Mighty  Epic  Model) 


Shell: 


Steel  E  =  50) 
n 

a  =  lm  (39.3701  inches) 


Excitation: 


Internal  Triangular  Pulse 


Internal:  p(9,t)  =  Pq  f(t) 


t  — U 
r  f 


t^  is  a  convenient 
delay  constant 


Large  Yield 


Profile  1:  t  /t, 
r  f 


2,  t£/tA  =  8 


Near  Surface 
Small  Yield 


Profile  2:  t  /t, 
r  S 


0.2,  tf/t^  =  1.8 


Two  pressure  magnitudes:  PQ/Pe  =  PD/Pe  *  5 

where  Pe  is  the  maximum  pressure-loading  for 
which  the  medium  remains  elastic. 


Thus  we  have: 


Four  Internal-Excitation  Axisymmetric  Cases 


Table  2.  Material  Properties 


Inelastic  Medium  NTS  Tull  (Mighty  Epic  Model) 


Elastic  Portion 

K  =  100  Kb  =  1,470  ksi 
G  =  38  Kb  *  559  ksi 

Failure  Surface 

A  =  0.5  Kb  -  7.36  ksi 
B  =  0.52  Kb"1  =  0.03537  ksi"1 
C  =  0.44  Kb  =  6.486  ksi 

TOUT  =  0.1  Kb  =  1.47  ksi  (tension  cutoff) 


Cap 

R  =  3 

W  =  0.015 

X  =  -2.676  ksi 
o 

Mass  Density 

p  =  2  gm/cm3  =  3.879  slug/ft3  (125  ///ft3) 

Wave  Speeds 

C  =  9,100  fps 
P 

C  =  4,560  fps 
s 

Elastic  Steel 


E  =  30,000  ksi 
G  =  11,538  ksi 
v  =  0.3 

p  =  14.3  slug/ft3  (460  ///ft3) 

C  =  30,166  fps 
P 

Cg  =  10,780  fps 
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Table  3.  DAA-SSI  Check  Runs:  2 


Medium: 

Shell: 

Excitation: 


NTS  Tuff  (Mighty  Epic  Model) 

Steel  (a/h  =  50,  where  a  =  1  meter  =  3.28  feet) 

Internal  Triangular  Pulse,  Axisymmetric  Case 

Internal  Pressure  P(t)  =  P  f(t) 

o 

where  Pq  =  -jKb  =  7.35  ksi  and  f(t)  is  defined  as 
follows: 


Table  4.  Material  Properties 


Inelastic  Medium  NTS  Tuff  (Mighty  Epic  Model) 


Elastic  Portion 

K  =  100  Kb  =  1,470  ksi 
G  =  38  Kb  =  559  ksi 

Failure  Surface 

A  =  0.5  Kb  =  7.36  ksi 

B  =  0.52  Kb-1  =  0.03537  ksi-1 

C  =  0.44  Kb  =  6.468  ksi 

D  =  1.8  Kb-1  *  0.12245  ksi-1 

TOUT  =  0.1  Kb  =  1.47  ksi  (tension  cutoff) 


Cap 

R  =  3 

W  =  0.015 

X  =  -0.003  Kb  =  -0.441  ksi 
o 

Mass  Density 

Pm  =  2  gm/cm3  =  3.879  slug/ft3  (125  #/ft3) 

Wave  Speeds 

C  =  9,100  fps 
pm  ’  r 

C  =  4,560  fps 
sm  ’  * 

Elastic  Steel 


E  =  30,000  ksi 
G  =  11,538  ksi 
v  =  0.3 

pg  =  14.3  slug/ft3  (460  #/f t3) 

C  =  20,166  fps 
ps 

C  =  10,780  fps 
ss 
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SECTION  IV 


CONCLUSION 


The  numerical  results  of  the  previous  section  indicate  that  the  doubly  asymp¬ 
totic  approximation  produces  reasonable  results  for  nonlinear  medium-structure 
interaction  when  the  medium  obeys  ?  work- hardening  ^-plasticity  theory.  For  a  more 
realistic  soil  model  (viz.,  the  Cap  model),  however,  reasonable  results  have  not 
been  obtained.  In  particular,  the  aspect  of  the  Cap  model  that  creates  difficulties 
for  the  quasi-static  DAA  is  the  hydrostatic-stress  invariant  Ij.  This  impediment 
led  to  the  introduction  of  the  volume-wave  model  to  create  the  quasi-dynamic  DAA. 
Unfortunately,  the  dilatancy  and  tensile  cutoff  aspects  of  the  Cap  model  create  dif¬ 
ficulties  that  are  not  resolved  by  the  quasi-dynamic  DAA.  An  alternative  to  the  DAA, 
an  acceleration  body-force  approximation  based  on  the  volume-wave  model,  was  devel¬ 
oped.  This  approximation,  however,  possesses  marginal  numerical  stability  and  yields 
results  that  are  no  better  than  those  produced  by  the  quasi-dynamic  DAA. 

Perhaps  the  primary  conclusion  that  may  be  drawn  from  this  study  is  that  the 
uniformity  assumptions  inherent  in  boundary-element  formulations  cannot  be  avoided. 
Hence,  the  appearance  of  dilatancy  or  tension  cutoff,  which  lead  to  highly  nonuniform 
material  behavior,  cannot  be  tolerated  outside  the  surface  to  which  the  DAA  is  ap¬ 
plied.  This  implies  that  a  DAA  may  be  used  outside  of  a  nonuniform  region  treated 
with  finite-difference  or  finite-element  methods.  If  hydrostatic-stress  effects  are 
negligible  outside  that  region,  a  quasi-static  DAA  may  suffice.  If  they  are  not 
negligible,  a  volume-wave  model  may  be  required  to  produce  a  quasi-dynamic  DAA.  In 
the  event  that  internal-forcing  problems  are  satisfactorily  treated  in  this  manner, 
problems  involving  excitation  by  incident  waves  require  the  use  of  a  wave  decoupling 
approximation  such  as  that  discussed  in  Appendix  A. 

It  is  clear  from  the  preceding  discussion  that  boundary-element  and  infinite- 
element  methods  for  nonlinear  soil-structure  interaction  analysis  are  in  an  embry¬ 
onic  state  of  development.  Furthermore,  because  of  the  complexities  involved  in 
their  formulation  and  implementation,  the  prospects  for  near-term  utilization  are 
not  very  promising.  However,  the  payoff  associated  with  their  successful  applica¬ 
tion  is  sufficiently  great  that  promising  ideas  for  their  development  should  be 
pursued. 
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APPENDIX  A 


WAVE  DECOUPLING  APPROXIMATION 

This  appendix  presents  a  brief  discussion  of  the  approximation  introduced  into 
the  DAA  from  the  decomposition  of  surface  forces  and  displacements  into  incident  and 
scattered  wave  components.  A  simple  bar  model,  which  illustrates  the  physics,  is  con 
sidered. 

The  basic  concept  underlying  the  soil-structure  interaction  model  presented 
here  is  the  replacement  of  the  effectively  infinite  volume  of  medium  surrounding  the 
structure  by  an  interface  surface  that  provides  to  the  structure  an  approximate  repre 
sentation  of  the  dynamic  behavior  of  the  medium.  This  representation  is  one  that  ap¬ 
proaches  exactness  in  both  the  low-frequency  and  high-frequency  limits;  hence  the 
name  "doubly  asymptotic  approximation".  An  associated  concept  is  the  decomposition 
of  surface  forces  and  displacements  into  two  sets,  one  pertaining  to  the  incident 
wave  (fjjUj)  and  the  other  pertaining  to  the  scattered  wave  (fs>ug).  In  this  connec¬ 
tion,  a  second  approximation  is  introduced,  viz.,  a  "wave  decoupling  approximation" 
(WDA),  which  will  now  be  described. 

Figure  A-l  shows  a  rigid  mass  embedded  in  a  grounded  vertical  bar  that  is  load¬ 
ed  at  the  top  by  a  slowly  varying  "incident"  force  f  .  The  bar  material  has  zero  mass 
density  and  exhibits  bilinear-hysteretic  constitutive  behavior.  The  graph  labeled 
"exact"  shows  the  stress-strain  trajectory  for  a  point  just  below  the  mass  at  eight 
incremental  stages  of  the  loading  process.  The  graph  labeled  "approximate"  shows  the 
corresponding  trajectory  produced  by  the  WDA,  as  follows.  A  load  increment  Af^  is 
applied  to  the  bar  with  the  "scattered"  force  fg  =  Mg  neglected;  this  leads  to  Point 
1'  on  the  trajectory.  The  force  fg  is  then  considered,  which  leads  to  Point  1  on  the 
trajectory.  An  additional  load  increment  is  applied  to  the  bar  with  fg  neglected; 
this  leads  to  Point  2'.  Then  the  force  fg  is  again  considered,  which  leads  to  Point 
2.  This  process  continues  until  Points  6'  and  6  are  reached;  a  load  increment  Af^ 
is  then  subtracted  from  fj.  ?=  6Afj,  which  leads  to  Points  7'  and  7.  Finally,  another 
load  increment  Af^  is  subtracted  from  f^  =  5Afj,  which  leads  to  Points  8'  and  8. 

The  trajectories  of  Figure  A.l  indicate  that  the  WDA  is  in  error  only  during  un¬ 
loading,  where  path-dependent  constitutive  behavior  becomes  manifest.  This  implies 
that  the  WDA  should  be  accurate  during  the  loading  phase  of  a  multi-dimensional  dyna¬ 
mic  interaction  analysis  for  a  structure  whose  bulk  stiffness  and  mass  properties  are 
comparable  to  or  greater  than  those  of  the  surrounding  medium.  In  any  case,  the  WDA 
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sho-ld  be  accurate  uhen  f  is  much  snaller  thaa  f  .  because  f  . 

creasy  distance  fro-  the  surface  of  a  2-D  or  3  l  J  5  Ult'’  In- 

therefore  be  obtained  by  noving  the  interface  surface 

““  —sad  as  part  of  1  strucllre 
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APPENDIX  B 


OUTLINE  OF  THE  DERIVATION  OF  THE 
BOUNDARY  INTEGRAL  EQUATION  FOR  INELASTIC  MEDIUM 

This  appendix  outlines  the  initial  steps  in  the  derivation  of  the  boundary  in¬ 
tegral  equation  for  an  inelastic  medium  (5).  The  derivation  is  carried  out  to  the 
point  at  which  o^ther  investigators  [3,4]  have  incorrectly  neglected  a  term,  as  point¬ 
ed  out  by  Mukerjee  [5].  The  steps  that  involve  the  boundary  integral  equation  itself 
are  not  presented  as  the  use  of  Kelvin's  solution,  Betti's  reciprocal  theorem  and  the 
divergence  theorem  to  obtain  this  equation  has  been  thoroughly  studied  [3,4,5,22,23]. 

The  total  strain  e . .  is  written  as  the  sum  of  the  elastic  and  plastic  strain  to 
ij 

give 


e  p 

e  .  .  =  e  .  .  +  e  .  . 

13  13  13 


(B-l) 


The  total  strain  is  related  to  the  displacement  by  the  linear  kinematic  relation 

(B-2) 


ij 


<ui,j  +  uj,i>/2 


In  [3,4]  the  volumetric  plastic  strain  ej^  is  assumed  to  be  Eero;  an  incorrect  assump¬ 
tion  for  two-dimensional  plane  strain  [5]  and  for  three-dimensional  plasticity  theory 
involving  volumetric  changes,  such  as  the  Cap  model  [7],  Here  the  assumption 


ekk  *  0 


(B-3) 


is  made.  Hence,  from  Hooke's  law,  (B-l),  (B-2)  and  (B-3)  the  stress  becomes 


°ij 


Auk,k°"ij  +  y(ui,j  +  uj,i>  -  Aekk6ij  -  2vzh 


(B-4) 


where  X  and  p  are  the  Lame  coefficients.  Finally,  from  the  equilibrium  equation 


U  ,  .  .  —  —  U  . 

i3 .3  i 

and  (B-4)  the  Navier  equation  for  inelastic  medium  becomes 

u.  - i -  ..  =  2eP  +  _2v_  £P 


b  . 

l 


(B-5) 


(B-6) 


i , j j  (l-2v)  k , ki  “ij,J  '  l-2y  "kk.i 

where  G  =  p  and  Poisson's  ratio,  v  =  A/2(X+p). 

The  two  basic  paths  to  follow  in  the  boundary  integral  equation  derivation  start 
with  Betti's  reciprocal  theorem  and  use  (B-4)  [3]  or  start  with  B-6  [4],  In  both 
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V 


cases  the  underlined  term  is  not  included  in  [3,4],  but  Mukherjee  [5]  presents  a  cor¬ 
rect  derivation  of  the  kernels. 


As  the  only  difference  occurs  in  the  plastic  strain  term  the  kernels  for  the 
tractions,  displacements  and  body  force  in  (7)  remain  unchanged  from  the  linear  the¬ 
ory.  For  two-dimensional  plane  strain  they  are  [1,22,23] 


TkX 

T .  .  =  — 

ij  r 


ij 


nj 


akl  C4  +  2r 


i  j  .  k  r  i  j  ,  1  ) 


„  /  k  X  \ 

+  C.|n.  r..  .  -  n.  r..  , 

4 \  J  i J  ,  X  j  ij,k] 


(B-7  ) 


Xn  r  .  . 
ij 


i  j  ,  k 


1  J 
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and  the  correct  l  kernel  (8)  is  given  by  Mukerjee  [5]  as 


£  Xmk 
ij 
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mk 
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lj.k  ij  ,  X  ij,m 
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ij  .k 


(B-9) 


where  C^,  C^,  and  are  material  constants,  r_  is  the  distance  from  a  node  point 
on  the  i-th  element  to  the  variable  (field)  point  of  integration  on  the  j-th  element, 

k 

n^  is  the  unit  normal  to  the  surface  of  the  j-th  element,  n^  is  the  cosine  of  the 
angle  between  n^  and  the  k-th  Cartesian  direction,  6^  is  a  Kronecker  delta,  and  a 
subscript  following  a  comma  represents  spatial  differentiation  with  respect  to  the 
indicated  cartesian  coordinate  at  point  j . 


Equation  (B-9)  is  valid  for  two-dimensional  plane  strain.  Also,  the  £  kernel 
given  in  [3]  for  the  three-dimensional  case  is  incorrect  for  plasticity  which  in¬ 
cludes  volumetric  strain. 
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APPENDIX  C 


NUMERICAL  EVALUATION  OF  NONLINEAR  COEFFICIENT  MATRICES 


This  appendix  discusses  the  numerical  approach  used  to  evaluate  the  integral 

Jinik 

in  (8)  to  determine  the  matrix  elements  ,  The  evaluation  of  the  first  two  inte¬ 
grals  in  (7)  was  presented  in  [1]  and  the  integrals  in  the  third  of  (7)  and  (13)  and 
(14)  follow  the  approach  presented  here,  so  they  are  not  explicitly  considered. 

S-in 

As  the  shape  function  r\A  in  (8)  is  chosen  to  be  1  the  integral  (8)  reduces  to 

d  A  j  (C-l) 

j 

and  from  Appendix  B 
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where 


"ij  =  [(Xli  -  X13)2  +  U2i  '  X2j)Z]^ 


ij  ,k 


(xkj  -  xki)/rij> 


6^  j,  etc,,  are  Kronecker  deltas  and  =  1,2.  The  geometry  associated  with 

the  integration  is  shown  in  Figure  C-l.  This  illustrates  that  the  index  i  pertains 

to  a  fixed  node  point  and  the  index  j  pertains  to  a  point  within  the  quadrature  element 

Aj .  Note,  r^  is  the  distance  between  the  points  i  and  j  and  for  a  specific  r„, 

r..  ,  is  the  cosine  (k=l)  and  sine  (k=2)  of  the  angle  between  the  x, -global  axis  and 
x  j  » x 

the  line  r^.  The  quadrature  element  is  defined  as  an  eight  node  (parabolic  interpo¬ 
lation)  curvilinear  quadrilateral.  This  representation  allows,  where  applicable, 
direct  usage  of  finite  element  numerical  integration  methods  and  in  other  situations 
a  convenient  method  for  describing  the  geometry  through  the  shape  functions.  In  this 
example  the  plastic  strains  are  evaluated  at  the  center  of  element  (the  point  0  in 
Figure  C-l)  and  assumed  to  be  uniform  over  A^ ,  as  discussed  in  Section  2.2. 

Three  distinct  cases  arise  in  the  integration  of  (C-l):  1)  i  is  not  in  A^ 
(including  the  boundary)  and  a  ^  8  (see  Figure  C-l);  2)  i  is  not  in  A^  and  a  =  6  and 


C-l 


3)  i  is  in  A^ .  Case  1  is  easily  integrated  by  standard  numerical  integration  formu¬ 
las  over  the  square  shown  in  the  upper-righthand  corner  of  Figure  C-l.  That  is,  for 
integration  purposes,  the  curvilinear  quadrilateral  is  mapped  onto  the  square,  a 
standard  method  in  finite  element  analysis;  see  [24]  for  example.  Here  both  one- 
point  (C=n  =  0)  and  four-point  (C  =  n  discrete  integrations  were  found 

to  produce  3  to  4  figure  accuracy.  The  four-point  integration  rule  was  used  during 
the  latter  stages  of  the  study,  but  generally  the  one-point  integration  rule  is  suf¬ 
ficient.  For  r  (j  at  the  centroid  0)  greater  than  the  minimum  arc  length  from  i+1 
to  i  or  i  to  i-1,  the  general  rule  would  be  to  use  one-point  integration,  otherwise 
use  four-point  integration. 


Cases  2  and  3  are  covered  by  the  same  method.  Here  the  basic  problem  is  the 
singularity  in  the  integrand  (C-2)  for  case  3.  The  problem  for  case  2  appears  to 
come  from  the  even-odd  behavior  of  the  integrand  when  a  =  g.  The  method  that  re¬ 
solves  the  case  3  singularity  problem  also  resolves  the  case  2  difficulty,  so  they 
are  treated  together.  To  eliminate  the  singularity,  1/r^ ,  the  integration  over 
is  changed  to  polar  coordinates  p,8  with  the  origin  at  i  [25]  (see  Figure  C-2).  The 
r , 


...  ,  terms  are  functions  of  0;  see  discussion  below  (C-2),  so  (C-l)  now  has  the  form 

lj  ,k 

■Hi 


B  = 


f  (0)pdpd0 


(C-3) 


'e-'p 


Note  that,  with  this  polar  coordinate  system  p  =  r,  hence  (C-3)  becomes 


B  = 


SJ. 


f(8)dpd0. 


(C-4) 


so  the  singularity  is  removed.  The  integration  (C-4)  is  carried  out  over  the  four 
triangles  shown  in  Figure  C-2.  For  i  on  the  boundary  of  A^  only  3  or  2  triangles 
are  needed  to  cover  A^ ,  Also,  note  that,  f(8)  =  r^  lsee  equation  (C-2)].  The 

integral  (C-4)  is  evaluated  by  applying  the  three-point  one -dimensional  Gaussian 
quadrature  rule  [26]  first  to  dp  and  then  d8.  The  three-point  rule  is  used  as  it 
produces  accuracy  comparable  to  case  1. 
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integration  element  definition. 
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APPENDIX  D 


BOUNDARY  INTEGRAL  EQUATION  MATRICES 
WITH  SYMMETRY  CONDITIONS 


In  this  appendix  the  formulation  of  boundary  integral  equation  matrices  for  a 
symmetry  condition  is  presented.  To  illustrate  the  formulation  the  derivation  of 
the  linear  stiffness  matrix  for  a  symmetry  condition  is  considered.  The  other  ma¬ 
trices,  £  and  .E,  are  treated  similarly,  so  they  are  not  considered  here. 

The  formulation  is  best  illustrated  by  considering  the  simple  example  shown  in 
Figure  D-l,  the  X2~axis  being  an  axis  of  symmetry.  Matrices  are  to  be  formulated  for 
only  one-half  the  problem  such  that  the  symmetry  condition  is  properly  represented. 
The  node  points  are  selected  as  shown,  such  that  two  nodes  are  on  the  symmetry  axis, 
e.g.,  nodes  1  and  11  in  this  example.  Note,  that  to  obtain  the  same  response  from 
the  half  (symmetry)  model  as  from  the  full  model,  the  half  model  tractions  must  be 
one-half  the  full  model  tractions  but  the  displacements  are  equal  at  nodes  1  and  11. 
Thus,  the  sum  of  the  two  half  models  produces  the  full  model.  For  the  two  half  mod¬ 
els  the  matrix  equation  for  displacements  and  tractions  may  be  written  as 


where  the  subscript  r  indicates  the  retained  degrees  of  freedom  and  the  subscript  e 
indicates  the  degrees  of  freedom  eliminated  by  applying  the  symmetry  condition.  The 
submatrices  in  (D-l)  are  composed  of  the  following  elements: 

~11:  22  rows  and  columns  of  £  (Note,  there  are  22  rows  because 

there  are  11  nodes  and  2  DOF/node) 

S10:  The  first  22  rows  and  the  last  20  columns  of  S  and  the  last  2  columns 
~12  ~ 

of  jS^2  are  the  first  2  columns  and  22  rows  of  £. 

:  The  first  20  rows  and  22  columns  are  the  last  20  rows  and  first  22 

~21 

columns  of  S  and  the  last  2  rows  and  22  columns  are  the  first  2  rows 
and  22  columns  of  S. 

Of 

Sgj1  The  first  20  rows  and  columns  are  the  last  20  rows  and  columns  of  S, 
the  last  2  rows  and  20  columns  are  the  first  2  rows  and  last  20  col¬ 
umns  of  S,  the  last  2  columns  and  20  rows  are  the  first  2  columns  and 

iv  7 


the  last  20  rows  of  jS,  and  the  last  2  rows  and  2  columns  are  the  first  2 
rows  and  columns  of  jS. 

For  £i2’  ~21  aTl<*  ~22  t*le  arran8ement  of  elements  based  on  F  are  identical  to 

~11’  ~12’  ~21  an<*  -§22  wit^  an  exception.  The  exception  is  that  the  elements  of  the 
first  and  last  two  columns  of  each  of  the  four  submatrices  are  divided  by  two;  this 
accounts  for  the  one-half  value  of  the  traction.  In  the  above  _S  and  £  refer  to  the 
full  problem  matrices;  see  equation  (6). 

The  e-degrees  of  freedom  are  related  to  the  r-degrees  of  freedom  by 

u  =  T  u  and  t  =  T  t  (D-2) 

-e  ~  -r  -e  ~  -r 


where  the  transformation  matrix  T  imposes  the  symmetry  condition.  For  this  problem 
the  components  in  the  x^-direction  are  equal  in  magnitude  and  opposite  in  sign  and 
the  components  in  the  X2~direction  are  equal  in  magnitude  and  sign.  Thus,  T  is  a 
very  sparse  matrix  composed  of  +  and  -  l's. 

Enforcing  the  equivalence  of  virtual  work  for  the  full  model  and  the  symmetry 
model  gives 
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(D-3) 


where  superscript  T  indicates  matrix  transposition.  Equations  (D-3)  may  appear  to 
require  as  much  storage  for  assembly  as  the  full  model  matrices,  but  they  do  not  as 
each  of  the  four  terms  are  assembled  individually  and  stored  in  and  as  they 
are  generated.  The  stiffness  matrix  for  the  symmetry  model  is  then  formed  from 

K  =  F-1  S  (D-A) 


If  a  method  existed  to  directly  generate  the  boundary  element  stiffness  with¬ 
out  forming  F  and  £  and  then  solving  for  K  the  matrix  algebra  presented  could  be 
eliminated.  The  K  matrix  would  simply  be  assembled  for  half  the  problem  and  the 
boundary  conditions  for  symmetry  imposed,  as  is  done  in  the  finite  element  method. 
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Rome  Air  Development  Center 
Air  Force  Systems  Command 
ATTN:  TSLD 

Strategic  Air  Conmand 
Department  of  the  Air  Force 

ATTN:  NRI-STINFO,  Library 

United  States  Air  Force  Academy 
ATTN:  DFCEM,  W.  Fluhr 


U.S.  Nuclear  Regulatory  Commission 

ATTN:  R.  Whipp  for  Div.  of  Security 
for  L.  Shao 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 


Aerospace  Corp. 

ATTN:  L.  Selzer 

2  cy  ATTN:  Technical  Information  Services 

Agbabian  Associates 
ATTN:  C.  Bagge 
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Art6C  Amci%tes- 

ATTN-  s.  Gin 

Avco  Research  g  cuc, 

ATTN-  iih^yStems  Group 
"•  Llbrary,  fl830  ^ 

BDM  Cor p. 

ATTN:  Corpora  te  Library 
BOM  Corp, 

ATTfi:  B.  Hensley 

Be"  ATT?*0"*  Lat>S 
AT™:  J.  White 

Boeing  Co. 

Iff  f’S&V"’"0’ 

TN-  J.  Wooster 

California.  Insti tute^of  Technology 
California  Research  *  r  . 

Sff  «•  4i£20h'®’-  "*• 

^"ras**..  <*. 

'"'"asr-ife*. 

gteiisjr 

zee,  Officer  for  i  u- 

EMr  u  i..  d<  WlsPtski 

“r's; *«.«.  «**, 

flTTN:  s.  SI  iter 
electromechanical  Svs  u 

ATTN:  «•  Shook  f  NSW  the. 

-fric,».  Wang 

*  Ca)houn 
franklin  institute 
AT™:  z.  Zuda„s 

^ra?  Dynamics  Corp. 

ATTN:  K.  Anderson 


,c 

- LContu^ 


Genera)  nectric  Co. 

ATTN:  M.  Bortner 

Genenlnlectr <c  Co. 

ATTN:  A.  Ross 
Genera  1  Electric  r„ 

ACTN:  MsjAC  y'TEMP° 
General  Research  Corp_ 

N-  B.  Alexander 
Geocenters,  inc 

AT™:  E.  'Marram 
H-Tech  tabs,  Inc 

AJTN:  B.  Hartenbaum 

111  RTmCh  r,nsti  tute 

ATTN*  .0cumen*$  Library 

w*  Longinow 

H  far* 

0.  HaltTwanger 

""Avar 

J'  JS®"*,0’;.  «*. 

A1TN:  J.  Collins 
liaman  Avi  Dyne 

ATTN:  Library 

l<ama"  Ixlfnces  c°rP. 

ATTN:  Library 

ICara90A"nNn.and  Case 

•  J.  Karagozi an 

Lockheed  Missiles  g  rna  „ 

ATTN:  ZLC,  Library  ;"c- 

Lockheed  Missiles  and  s 

ml:  ?•  Co-  '*■ 

AUh  J'  Geers 
•  p.  Underwood 

C°rP. 

ATTN:  A.  Cowan 

S3*. 

McDonn,e"  Gouglas  Corp_ 

AT™:  R.  Halprin 
^Titt  CASES,  lnc. 

Si  Library* 
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DEPARTMENT  OF  DEFENSE  CONTRACTORS  (Continued) 


Meteorology 

Research,  Inc. 

ATTN 

W.  Green 

Mitre  Corp 

ATTN 

Di rector 

Uni  verst  ty 

of  Oklahoma 

ATTN 

J.  Thompson 

Pacifica  Technology 

ATTN 

G.  Kent 

ATTN 

R.  Bjork 

ATTN 

R.  Allen 

Physics  International  Co. 

ATTN 

E.  Sauer 

ATTN 

L.  Behrmann 

ATTN 

£.  Moore 

ATTN 

Technical  Library 

R&D  Associates 

ATTN:  R.  Port 
ATTN:  C.  MacDonald 
ATTN:  J.  Lewis 

ATTN:  Technical  Information  Center 

Rand  Corp. 

ATTN:  Library 
ATTN:  C.  Mow 

Science  Applications,  Inc. 

ATTN:  Technical  Library 

Science  Applications,  Inc. 

ATTN:  S.  Oston 

Science  Applications,  Inc. 

ATTN:  R.  Hoffmann 
ATTN:  D.  Bernstein 

Science  Applications,  Inc. 

ATTN:  B.  Chambers,  III 

Southwest  Research  Institute 
ATTN:  W.  Baker 
ATTN:  A.  Wenzel 


DEPARTMENT 

OF  DEFENSE  CONTRACTORS 

SRI  International 

ATTN 

G.  Abrahamson 

ATTN 

W.  Wilkinson 

Systems,  Science  X  Software,  Inc. 

ATI  N 

Library 

ATTN 

T.  Cherry 

ATTN 

T.  Riney 

ATTN 

D.  Grine 

Teledyne  Brown  Engineering 

ATTN 

J.  Ravenscraft 

Terra  Tek, 

Inc. 

ATTN 

Library 

Tetra  Tech 

Inc. 

ATTN 

Library 

Texas  A  X  M  University  System 

Texas  A  X  M  Research  Foundation 

ATTN 

H.  Coyle 

TRW  Defense  X  Space  Sys.  Group 

ATTN 

Technical  Information 

ATTN 

P.  Bhutta 

ATTN 

A.  Narevsky 

2  cy  ATTN 

P.  Dai 

Center 


TRW  Defense  &  Space  Sys.  Group 
ATTN:  G.  Hulcher 


Weidlinger  Assoc.,  Consulting  Engineers 
ATTN:  J.  McCormick 
ATTN:  M.  Baron 

Weidlinger  Assoc.,  Consulting  Engineers 
ATTN:  J.  Isenberg 

Westinghouse  Electric  Corp. 

ATTN:  W.  Volz 
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